Cluster Synchrony in Systems of Coupled Phase Oscillators with Higher-Order Coupling 
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We study the phenomenon of cluster synchrony that occurs in ensembles of coupled phase oscillators when 
higher-order modes dominate the coupling between oscillators. For the first time, we develop a complete analytic 
description of the dynamics in the limit of a large number of oscillators and use it to quantify the degree of 
cluster synchrony, cluster asymmetry, and switching. We use a variation of the recent dimensionality-reduction 
technique of Ott and Antonsen (Chaos 18, 037113 (2008)) and find an analytic description of the degree of 
cluster synchrony valid on a globally attracting manifold. Shaped by this manifold, there is an infinite family 
of steady-state distributions of oscillators, resulting in a high degree of multi-stability in the cluster asymmetry. 
We also show how through external forcing the degree of asymmetry can be controlled, and suggest that systems 
displaying cluster synchrony can be used to encode and store data. 
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I. INTRODUCTION 



Large systems of coupled oscillators occur in many exam- 
ples throughout science and nature and serve as a basic model 
for emergent collective behavior. Examples include synchro- 
nized flashing of fireflies (TJ], cardiac pacemaker cells y|], 
walker-induced oscillations of the Millennium Bridge yD, and 
circadian rhythms in mammals J4}]. Under certain conditions, 
these limit cycle oscillators can be approximately described 
entirely in terms of their phase angles 9. Kuramoto showed ||5[] 
that the evolution of the phases in an ensemble of TV weakly 
coupled oscillators obeys 
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where 6 n and uj n are the phase and intrinsic frequency of os- 
cillator 71, and H nm is a 27r-periodic function. The choice 
of H nm (0) = (K/N) sin(#), which leads to the Kuramoto 
model |5J], is well motivated because the expansion of sev- 
eral coupled oscillators about a Hopf bifurcation generically 
leads to sinusoidal coupling. This choice has also become 
a paradigm for the study of emergence of synchrony in cou- 
pled heterogeneous oscillators |gj. Many generalizations of 
the Kuramoto model have been studied. Some examples are 
systems that account for time-delay ]7[], network structure Ja], 
non-local coupling J3], external forcing IUOI]. non-sinusoidal 
coupling lUlll . Josephson junction circuits IU2I1 , coupled ex- 
citable oscillators IU3I1 . bimodal distributions of oscillator fre- 
quencies IU4I1 . phase resetting IU5I1 . time-dependent connec- 
tivity JUa], noise IU7I1 . and communities of coupled oscilla- 
tors IU8I1 . Recent analytical work IU9H21TJ (in particular the 
Ott-Antonsen (OA) ansatz IU9I0 has allowed for the simplifi- 
cation of the analysis of these systems to the study of reduced 
low-dimensional equations and made many of these systems 



analytically tractable. 

While the choice H nm (0) = (K/N) sin(0) that yields the 
Kuramoto model is the simplest, describes many situations of 
interest, and has the advantage of being analytically tractable, 
more general choices can result in additional dynamical fea- 
tures. If there are higher harmonics in H nm , but the sinusoidal 
term is dominant, there is a transition to synchrony as in the 
Kuramoto model as the coupling between the oscillators in- 
creases llllll . In this case, the synchronous state is character- 
ized by the phases of a large fraction of the oscillators clus- 
tering around a common phase. When higher harmonic terms 
are dominant, however, the synchronous state is characterized 
by the formation of multiple synchronized groups (or "clus- 
ters") of oscillators, each with a common phase 122I1 . This 
phenomenon has also been called multibranch entrainment in 
previous work 112311 . Cluster synchrony occurs in many ap- 
plications in nature, including networks of neuronal, photo- 
chemical, and electrochemical oscillators II24h26|1 . as well as 
genetic networks J27I1 . In this paper we will study Eq. ((TJ with 
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for integer q > 2, which, as we will see, leads generi- 
cally to the formation of q clusters. There are various ex- 
perimental and theoretical motivations for the study of this 
model. In Ref. (25( 1. experiments with systems of globally 
coupled photochemical oscillators were performed in which 
two synchronized clusters emerged. In Ref. 112611 . the cou- 
pling function between electrochemical oscillators n and m 
was directly measured and found to be qualitatively equiv- 
alent to either H nm (6) = (K/N)sm{9) at a lower volt- 
age, which is equivalent to the classical Kuramoto model, or 
H nm (6) = (K/N) sin(20) at higher voltage, which is equiv- 
alent to Eq. (O with q = 2. In some Kuramoto-type models of 
neuronal networks, a coupling function of the form in Eq. (O 
and the associated cluster synchrony arises as a result of learn- 
ing and network adaptation. It has been proposed that the cou- 
pling between oscillators in such networks evolves according 
to a Hebbian learning rule 1280. If this learning is fast, Eq. (O 
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FIG. 1: (Color online) Example oscillator configurations of a system with (a) no synchrony (|ri|, |?"2 1 ~ 0), (b) symmetric (|n| ~ 0) cluster 
synchrony (|t"2 | > 0), and (c) asymmetric (|n| > 0) cluster synchrony (|r-2 1 > 0) and corresponding density functions (d), (e), and (f). 

fraction of oscillators switching to a cluster around 6 = $. If 
different values of 'J 7 are chosen for different oscillators (i.e., 
"J (->• $> n ), then if F is large enough with respect to \u> n \ and 
K, the phase of oscillator n will converge to a value near ^f n . 
This paper is organized as follows. In Section II we solve 
for the dynamics of the system with Eq. (O and q = 2. In Sec- 
tion III we discuss the effect of external forcing on asymmet- 
ric clustering and switching, the presence of hysteresis when 
the coupling strength is changed, as well as how asymmetric 
clustering can be used to store information. In Section IV we 
discuss how results generalize to the case q > 2. Finally, in 
Section V we conclude this paper by discussing our results. 



is recovered with q = 2 [29]. We note that the applications 
mentioned above all use q = 2, but larger q values are also 
relevant. For instance, cases of three or more clusters have 
appeared in the study of synthetic gene networks 12711 . 

Cluster synchrony has been studied in many contexts, for 
example in networks of phase oscillators with 112211 and with- 
out noise OCX 13111 . networks of integrate-and-fire oscillators 
113211 . and more general cases 0311 . Reference J28ll studied 
Eqs. (Q3 and (O in steady-state using a self-consistent ap- 
proach to characterize the phase distribution and stability of 
the clusters. Reference 112211 studied the dynamics of clusters 
in ensembles of oscillators when the coupling function H has 
two Fourier modes under the effect of small noise. Despite 
these and other studies IBOIl . a complete analytical treatment of 
Eqs. (Q]i and ® is lacking. For example, Ref. [28] studies the 
steady state solution using a self-consistent approach, but does 
not analyze the dynamics, while Refs. 1221 13011 assume iden- 
tical oscillators. In this paper we will use the Ott-Antonsen 
ansatz to obtain a low-dimensional description of cluster dy- 
namics and a full solution to Eqs. (Q]) and (0. Thus, our so- 
lution of Eqs. (fl} and (O is analogous to the recent solution 
H19L I2OII of the Kuramoto model in that, even though partial 
solutions existed previously, our solution fully characterizes 
the dynamics (with the same caveats as in Refs. IU9ll20lD . 

Two interesting phenomena that are particular to systems 
displaying cluster synchrony are asymmetric clustering [13411 
and switching 1221 12511 . Asymmetric clustering is character- 
ized by a non-uniform distribution of oscillators in different 
clusters and switching refers to oscillators moving between 
clusters. We find that asymmetric clustering emerges from 
non-uniform initial conditions, to which systems with a cou- 
pling function of the form of Eq. (|2]) with q > 2 are very sen- 
sitive. Switching can be achieved by introducing an external 
forcing term, Fsin(\I/ — cu t — 8 n ) (where cj is the average 
oscillator frequency), on the right hand side of Eq. (Q]) with 
F y^ nonzero for a finite amount of time. This results in a 



II. LOW-DIMENSIONAL DESCRIPTION OF THE 
TWO-CLUSTER STATE 

In this section, we will study in detail Eqs. (Q]) and (0 with 
q = 2, which leads to the system 
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where the intrinsic frequencies u n are drawn randomly from a 
distribution g(co). Also, we define the set of generalized order 
parameters 



a- 



r k e 



i^k 



1 

N 



N 
m— 1 



Akdr, 



(4) 



for k £ N. These generalized order parameters were intro- 
duced in previous work 113511 where coupling functions with 
higher harmonics were studied. We will see that when more 
than one cluster emerges, more than one rk is needed to de- 
scribe the dynamics of the system. In this case of q — 2, two 
clusters emerge. The order parameter |ra| measures the de- 



gree of cluster synchrony in the system while \r±\ measures 
the degree of asymmetry in clustering (see Fig.[TJ. Eq. can 
be rewritten in terms of r-i as 
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where * denotes complex conjugate. 



Following Ref. I19H . we let N — > oo and move to a con- 
tinuum description. Accordingly, we introduce the density 
function f(9, u, t), which describes the density of oscillators 
with phase 9 and natural frequency lo at time t. Since os- 
cillators are conserved / must satisfy the continuity equation 
d t f + d e {fd) = 0, giving 



d t f + de 
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= 0. (6) 



To analyze Eq. @, we find it convenient to define the sym- 
metric and antisymmetric parts of /, f s and f a , as 

f 8/a (9,uj,t) = [f(9,uj,t)±f(9 + n,LJ,t)}/2, (7) 

where f s and f a are symmetric and antisymmetric with re- 
spect to translation by n, respectively, in the sense that f s (6 + 
ir,u,t) = f a (6,u,t)wdfa(0 + ir,w,t)=-fa(P,w,t). We 
note that / is a solution of Eq. (O if / = f s + f a and f s and f a 
are both solutions of Eq. ©, Thus, we can study separately 
the symmetric and antisymmetric dynamics of solutions /. 



A. Symmetric Dynamics 



We now use a variation of the OA ansatz to find a low- 
dimensional analytical solution for the dynamics of the sym- 
metric part of /, which evolves independently from the anti- 
symmetric part. For the Kuramoto model, after expanding the 
distribution / in Fourier Series, 
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where c.c. denotes complex conjugate terms, Ref. 111911 intro- 
duces the ansatz f n (u),t) = a n (ui,t) which yields a solu- 
tion for systems with sinusoidal coupling provided a evolves 
according to a simple ordinary differential equation (ODE). 
Solutions of this kind turn out to form a low-dimensional, 
globally-attracting, invariant manifold to which solutions con- 
verge quickly given that the spread in g(u>) is non-zero 112011 . 
This manifold is the set of Poisson kernels, 
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In terms of the Fourier series (|8), the symmetric part of / is 



given by 
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For the new system defined by Eq. (0, we use the follow- 
ing variation of the OA ansatz on the symmetric part of /: 
fan(w,t) = a n (u>,t). When Eq. ( fTOb with this ansatz is in- 
serted into Eq. (|6]l and projected onto the subspace spanned 
by e m6 , all equations reduce to the following ODE for a: 



a + 2iuia + K (i"2a — r^) = 0. 
In the continuum limit, we have 

/OC />27T 
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We now assume that g(cu) is Lorentzian with mean luq and 



spread A, i.e. g(u) 



Furthermore, by en- 



7T(A 2 + (W-W t 

tering the rotating frame 9 M> 9 + uiot we can assume without 
loss of generality that ujq = 0. With this choice of g{oS) we 
can integrate Eq. ( TT2b exactly by closing the contour with the 
semicircle of infinite radius in the lower-half complex plane 
and evaluating a at the enclosed pole (see Refs. lfl9l 12011 for 
the validity of this procedure): 



r s (t) = o*(-»A,t)=o*(t), 



(13) 



where we've defined a(t) = a(—iA,t) to simplify notation. 
By evaluating Eq. (fTTT i at cj = — iA, close the dynamics for 

T2- 



f 2 = -2Ar 2 + K[r 2 - r* 2 r\). 
In polar coordinates, r 2 = \r2\e 1 ^ 2 , we find 
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(15) 
(16) 



Thus, the unsynchronized state (i.e. [7^2 j = 0) is stable for 
K < 2A, at which point it loses stability and the stable syn- 
chronized branch |t"2 I = \/l — 2A/K emerges. 

We now show that solutions of the form given in Eq. ( fTOb 
with f n (w, t) = a n (ui, t), where a obeys Eq. (fTTT i. are globally 
attracting. An alternative way of solving for the dynamics of 
r2 is to make the change of variable <fi = 29, which yields a 
new continuity equation: 



dtf s +dd 
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which is of the same form of the equation studied in Ref. ||2C 
There it is shown that the dynamics of r 2 given by Eq. (fT4t 
are globally attracting provided that the spread of g(u>) is non- 
zero. Thus, the globally attracting invariant manifold for f s is 



the set of double Poisson kernels centered at i/j 2 , fs(0, u>, t) 
P (29 -^ 2 ,\r 2 {t)\,uj), where 



p(e, P ,u) 



9(w) 
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Since the system is invariant to rotations 6 i-> + y>, hereafter 
we assume without loss of generality that ip2 = 0. 



B. Steady-State Solution 



To find the steady-state solution /*f , we set dtf*J 



For |w| < K\r 2 \ we find that f*J a /g(co) = c hs/a S(6 - 
4>(oj)) + C2 lS / a S(6 — <I>(oj) — 7r), where <f>(oj) and 0(w) + it 
are the stable fixed points of Eq. (0 defined by c/)(uj) = 
i arcsin[w/(i ; tr|r2|)]. (Recall that we assume t/j 2 = 0.) Im- 
posing symmetric and antisymmetric conditions, we have that 
ci, s = c 2 , s = 1/2 and ci, a = -c 2 , a = c with |c| < 1/2. 



We first find the steady-state solutions of the system de- 
scribed by Eq. (0. Recall that the symmetric and antisymmet- 
ric parts of / satisfy the partial differential equation (PDE) 



d t f s/a + do [f s/a (uj ~ K\r 2 \sm(29))] . 
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For \lj\ > K\r 2 \, we find that / 
K\r 2 \sm(28)\, where C(w 



' /g{u>) = C(u> )/\u> - 
2^/uj 2 -K 2 \r 2 \ 2 /ir and 



fa S = 0- Thus, the complete steady-state distribution of os- 
cillators is 



ne,u) 



g(w)[(l /2 + c)6(0-<l >(w)) + (1/2 - c)5(9 - <f>(w) - it)] if |w| < K\r 2 \, 

2 9 (uj)^/lj 2 - K 2 \r 2 \ 2 /\ir[uj - K\r 2 \sm{29)}\ if \u\ > K\r 2 \, 



(20) 



r 



with I t 2 I = y/l — 2A/K. The interpretation of the different 
terms in Eq. ( f20T > is the following. For |w| < K\r 2 \ f ss is 
comprised of two delta functions representing the two clusters 
of phase-locked oscillators at 9 — <fi(uj) and (/>(w) + n. For 
|cj| > K\r 2 \ oscillators drift for all time and the second line 
in Eq. d20b is their steady-state distribution. 

While the symmetric part of the distribution is only depen- 
dent on the value of K, the antisymmetric part of the distri- 
bution depends on the free parameter c, which must be deter- 
mined from initial conditions. Thus, different solutions with 
different degrees of cluster asymmetry coexist. 

Now we compute the degree of cluster synchrony and 
asymmetry in the system at steady-state in terms of initial 
conditions. The degree of cluster synchrony is exactly \r 2 \ = 
y/l — 2A/K, but the degree of asymmetry, measured by | r\ | , 
depends on the free parameter c which must be determined 



from initial conditions. To calculate n, we note that only the 
locked portion (|w| < -ftT^D of / contributes to n, so 

/K\r 2 \ r2ir 
/ f ss {9,uj)e ie d9duj (21) 

-K\r 2 \ JO 
/K\r 2 \ 
ff (w)e l * ( "W (22) 

-K\r 2 \ 

Through a series of substitutions, this integral can be evalu- 
ated exactly: 



\n\ = 



2\/2c /arctan(A ) arctanh(A+)\ 

~V~ \ 1+ 1= ) ' (23) 



where A = 
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^K 2 \r 2 \ 2 + A 2 ±K\r 2 \ 
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FIG. 2: (Color online) Stable (solid blue) and unstable (dashed red) 
equilibria of 8 as a function of a; for phase-locked oscillators. Bound- 
aries between locked and drifting regions (w = ±A"jr2|) are plotted 
in black dotted lines. 



0.8 


lr 2 (t)l, b = 0, 0.4, 0.8 




0.6 


lr (t)l, b = 0.8 


E 0.4 


" lr (t)l, b = 0.4 




0.2 

n 


lr Ct>l, b = 



12 3 4 5 

time t 

FIG. 3: (Color online) Order parameters |fi(t)| (solid colored 
curves) and |f2(£)| (dashed colored curves) for 6 = 0, 0.4, and 0.8 
from a simulation of Eq. {3) with N = 10000 oscillators and ana- 
lytic predictions of steady-state (black dot-dashed lines). Parameters 
are K = 4, A = 1. 




FIG. 4: (Color online) Example characteristics 8(t) from Eq. dAlt 
for K — 4 and A = 1 of (a) locked oscillators (w = 1) and (b) 
drifting oscillators (w = 3). 

As an example illustrating the dependence of c on ini- 
tial conditions we assume for simplicity that the symmet- 
ric dynamics are at steady-state by time t = to so that 
|r*2 1 = \/l 2A/A", but the antisymmetric part may still 
not be at rest. Thus, phase-locked oscillators with natural 
frequency u n settle to one of the two stable equilibria of 
n = oj„. — K\r 2 \ sin(20 n ), while the unstable equilibria serve 
as boundaries for the basins of attraction. In Fig.|2]we plot the 
stable equilibria in blue solid lines and the unstable equilib- 
ria in red dashed lines for K = 4 and A = 1. Boundaries 
between locked and drifting regions, ui = ±A"|r2|, are plot- 
ted in dotted black lines. We denote the unstable equilibria by 
01,2 = —\ arcsin[w/(A'|r2|)] T f- From Eq. ( TSOl i we find 
that c + 1/2 is just the fraction of oscillators in the locked re- 
gion between 0i and 82, so c in terms of the initial density 
f(6, u), to) is 



To test this result, we choose initial conditions 
/(0,w,io)=P(20,p 2 ,w)[l + &cos(0)] ) 



(25) 



(26) 



which has symmetric and antisymmetric parts f s = 
P(20, |r 2 |,w) and /„ = 6P(26>, |r 2 |,w)cos(0), respectively. 

Choosing A' = 4 and A = 1 we integrate Eq. (f25t numer- 
ically to get c fa 0.4423516. In Fig. [3] we plot results from 
a direct numerical simulation of Eq. (0 compared with the 
analytical prediction above. We simulate N = 10000 oscilla- 
tors with K = 4 and A = 1 and plot |ri (£)| for b = 0, 0.4, 
and 0.8 in blue, red, and green solid lines (labeled in Fig. |3), 
respectively, with the predicted value of limt_}. 00 |ri(t)| for 
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FIG. 5: (Color online) Transient dynamics of |fi(i)| for initial con- 
ditions f(8, u),0) — P(29, \r2\,u))(l + bcos(8)) from simulation 
with N = 10000 and 6 = 0, 0.4, and 0.8 (blue, red and green 
curves) and from integrating Eq. J32t numerically (cyan, magenta, 
and black dashed curves). Parameters are K = 4 and A = 1. 

each in black dot-dashed. We also plot |?"2(t)| for each value 
and the predicted value of lim t _>. 00 |? , 2(t)| = l/v2 in black 
dashed curves. Simulations agree very well with the theory. 
Note that, unlike |ri |, |r2| (both predicted and from simula- 
tion) does not depend on b. 



C. Transient Dynamics 

From Fig.[3]we see that the \n \ dynamics reach steady-state 
quickly. To capture the transient dynamics we can solve the 
PDE© 



dtf 



K\r 2 \ sin(26)}def - 2A|r 2 | cos(20)/ (27) 



coupled with the \r 2 \ dynamics, which evolve independently, 
via the method of characteristics IB6I1 . The characteristic 
equations (along with ui = 0) are 



6» = w-A'|r2|sin(26'), 

/ = 2A|r 2 |cos(20), 
A, 
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(29) 

(30) 



When \r 2 \ is at steady-state Eqs. (f28lf30b can be solved an- 
alytically. Analytic expressions for the characteristic curves 
0(t,Oo) starting at the initial phase 0q and the distribu- 
tion f(6,ui,t), starting with initial condition f(0,cj,to) = 
g(uj)h(0) are given in Appendix A. 

Using K = 4 and A = 1, we plot some example charac- 
teristics using the analytic solution for u> = 1 and cj = 3 in 
Figs-Sta) andHfb), respectively. For these parameter values 
u = 1 is in the locked population and u = 3 is drifting. For 
u = 1, characteristics (solid colored lines) quickly converge 
to one of the two stable fixed points, with basins of attrac- 
tion separated by unstable fixed points (black dotted lines). 
Thus, / evaluated at uj = 1 converges very quickly to two 
point masses. However, for u = 3, the characteristics con- 
tinue drifting with a finite velocity for all time. 
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FIG. 6: (Color online) Numerically-computed distribution f(6,uj, t)/g(u}) obtained from numerically solving Eq. l |27b at times t = 0.33 (a), 
t = 0.67 (b), and t = 1 (c) with inital conditions P(20, 0.1, uj)(1 + O.4cos(0)) and parameters K = 4 and A = 1. 

In principle, we could calculate r\ (t) through the integral exactly [Qj 



n(t) 



f(6,uj,ty e d6duj, 



(31) 



— OO J — 7T 



where f(6,u>,t) is given by Eq. ( IA2I ) in Appendix A. How- 
ever, given the quick convergence of / to delta functions in 
the locked regime, Eq. OH is difficult to integrate numeri- 
cally, so we rather calculate r x (t) via the integral 



n(t) 



— OC J — 7T 



f{6{t,e Q ),u,ty e ^—dd du. (32) 
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In Fig. |5] we compare the results of integrating Eq. d32l l nu- 
merically with the simulations of Eq. (O using N = 10000 os- 
cillators, K = 4, A = 1, and/(0,w,O) = P(26, M,w)(l + 
6003(6*)). For 6 = 0, 0.4, and 0.8, |ri| obtained from simu- 
lations are plotted as solid lines, and results from integrating 
Eq. d32l l numerically are plotted as dashed lines. The results 
from Eq. d32l l capture the transient dynamics very well. 

The example above leading to Fig. [5]was for a case with |ra | 
initially at steady state. If |7~2 1 is not initially at steady-state, 
the solution to Eq. ( [TBI ) with initial condition |r"2 (0) | = po is 



hffll = P2I 
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PO 



a 2(2A-K)t 



(33) 



where P 2 = y/l - 2A/K. 

In Fig. [6]we plot the evolution of f(0, to, t) obtained from 
numerically solving Eq. d27l ) when the symmetric dynam- 
ics are not at steady-state. Starting with initial conditions 
/(0,w,O) = P(26, 0.1, w) (1 + 0.4 cos{0)) and parameters 
K = 4, A = 1 we plot the distribution/ (0, u),t)/g(u) at 
t = 0.33 (a), t = 0.67 (b), and t = 1 (c). We see that the dis- 
tribution quickly localizes, in agreement with the asymptotic 
form in Eq. d20l . In Fig.|7]we compare |ri(t)| and ^(t)! cal- 
culated from the numerical solution of Eq. d27b (blue circles 
and red triangles) with the same variables calculated from a 
direct simulation of Eq. ([3]) with A^ = 10000 oscillators (cyan 
and magenta dashed lines). The analytic solution |r*2 (*) | in 
Eq. ( f33l > is plotted as a black dot-dashed line. 



III. EXTERNAL DRIVING AND HYSTERESIS IN THE 
TWO-CLUSTER STATE 




time t 

FIG. 7: (Color online) Comparison of |n | and |ra | from the numeri- 
cally solved PDE Ml\ (blue circles and red triangles) and from direct 
simulation of Eq. (O with iV = 10000 oscillators (cyan and magenta 
dashed lines) with the same initial conditions and parameters used in 

Rg. m 



As we have seen, Eq. ([3]) admits a family of steady-state so- 
lutions characterized by a free parameter c. In this Section, we 
demonstrate that, by appropriately forcing Eq. OJ and modu- 
lating the coupling strength, the system can be driven to any of 
these solutions, thus allowing us to encode any desired value 
of c in the state of the system. Assuming luq = 0, we consider 
the forced system 



r^ N 



'n ^n 



N 



where F(t) 



Y, sin[2(0 m - 9„)} + F(t) sin($„ - B n ) 

rn— 1 

(34) 

F ifteJ 
otherwise, 



for some forcing magnitude Fq and time interval I = \t\, <a]- 



For Fq sufficiently large in comparison with \u n \ and K 
and duration d = £2 — ti not too small, 9 n will approach 
sa arcsin(w„/F ) + $„«$„ for F > |w„| + K. Thus, if 
$„ = for all n with d and Fo large enough (i.e. Fo 3> 
K\J\ — 2A/K), all locked oscillators are entrained to the 
= cluster and remain there after t = t<x, thus creating 
a completely asymmetric cluster state. On the other hand, if 
$„ are drawn from the distribution /i($) = d<5($) + (1 — 
d)6($> — 7r), then the ratio of the number of oscillators end- 
ing up in the cluster centered at ip = to those in the cluster 
centered at ^1 — n is d/(l — d), which forces c in Eq. (f23b to 
d— 1/2. Thus, by choosing appropriately the external forcing, 
we can set any degree of asymmetry we wish. 

To explore the effect of different Fq and d values, we sim- 
ulate Eq. (f34t with N = 2000 oscillators with random initial 
conditions and parameters K = 4 and A = 1 until steady- 
state (and attaining two clusters of approximately equal size, 
|t"x I ~ 0), then force the system with a strength of F for a 
duration d and all $ n = 0, then allow the system to reach 
steady-state and plot the resulting |?*i| value in Fig. |HJa). For 
very small Fo or d, \r-y\ remains small, but as soon as both are 
large enough the resulting \n | increases quickly. 

By forcing the system in this manner we achieve switching, 
i.e. oscillator n switches to the cluster centered at phase $„ 
if uj n is not too large. We note here that this kind of forced 
switching is qualitatively different than that in Ref. 112511 . In 
our original system given by Eq. (01, switching does not occur 
spontaneously. Thus, external forcing is necessary to observe 
the phenomenon. However, in Ref. J25I1 switching occurs 
spontaneously due to a heteroclinic orbit between different 
cluster states. 

Next, we consider the effects of slowly (compared with 
A -1 ) changing the coupling strength K after a steady state 
with some asymmetry is reached. If steady state is reached 
at t = to with a coupling strength K = Kq, then consider 
changing K to K\. We find hysteretic behavior in |ri| but not 
[7*2 1. Regardless of whether K\ < Kq or vice-versa, |ra| con- 
verges quickly to the predicted value [7*2! = \J\ — 2A/Ki, 
but the dynamics of |n| are more interesting: if K\ < Kq 
then 1 7"i I decreases significantly, but if K\ > Kq, then |n| 
remains approximately constant. In this situation, at time 
to, the distribution of oscillators is given by Eq. d20l i. If 
K\ < K the locked population loses all oscillators with 
yjKl - 2AF'i < |w| < yjKl - 2AK and \r x \ changes 
accordingly (maintaining the same c value, since these oscil- 
lators are lost in equal proportions from both clusters). On the 
other hand, if K\ > Kq the locked population will gain os- 
cillators with y/Ky - 2AK < \u\ < \J K\ - 2AK\. How- 
ever, at t = to the distribution for these drifting oscillators is 
perfectly symmetric, so both clusters pick up an equal number 
of oscillators and the symmetric density f a changes, while the 
antisymmetric density f a remains the same. Thus, the only 
change in \r\\ comes from the slight tightening of the phases 
4>{uS) = I arcsin[a;/(Ar|r i 2|)] and 4>{oS) + it about the clusters 
at 9 = and it. 

Extending this analysis to the case where K is both in- 
creased and decreased, |n| will never increase significantly, 
and only decrease significantly when K is decreased below a 
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FIG. 8: (Color online) (a) Steady-state |t*i | after forcing a symmetric 
distribution with forcing magnitude Fo for a duration d with K = 4 
and A = 1. (b) Hysteretic behavior of |n | (blue circles) vs the non- 
hysteretic behavior of \r2\ (black triangles) when K is changed in 
time (red dashed line). 

previous minimum. In figure [8{b) we plot |ri| and |ra| (blue 
circles and black triangles, respectively) as we change K (red 
dashed line). While |t~2 | follows the predicted behavior (green 
dot-dashed line) without any hysteresis, it is clear that \r\ \ be- 
haves as described above. 



We now suggest, as others have 0711 . that systems such as 
that given by Eq. (01 provide ways for encoding and storing 
data. These systems have the unique property that the sym- 
metric dynamics have a unique, (globally) stable fixed point, 
while there is a high degree of multi-stability in the antisym- 
metric dynamics. Furthermore, we have demonstrated above 
that through forced switching and modulation of the coupling 
strength, the asymmetry (i.e., |ri|) can be controlled. Thus, 
we suggest that a continuous valued variable could be stored 
and retrieved by representing it by |n|. Furthermore, in the 
general q case, which we study next, we will see that in addi- 
tion to one globally-attracting symmetric part, there are q — 1 
additional modes that display multi-stability. Thus, through 
similar techniques the q — 1 quantities |n(t)|, . . . , |r g _i(i)| 



can be controlled and used to store and retrieve q 
continuous valued variables. 



1 different 



IV. GENERAL g > 2 



We now discuss how the dynamics of the two state case 
generalize to higher-order coupling functions. Thus, we study 



the system 



K 

N 

K 



N 



^ sin[g(6» m - 9 n )] 



'r q e- qi9n - r*e qW " 



(36) 



for integer q > 2 and uj n randomly drawn from the distribu- 
tion g(uj). We find in this situation that q clusters form. 

Again, we introduce a continuum description and repre- 
sent the distribution of oscillators with the density function 
f(9, u, t), which satisfies the continuity equation 



dtf + Og 



f 


^ + |(r g e-^- 


-r,V»)) 



0. (37) 



In analogy with Eq. ( fTOb we define the modes 

9-1 

<1 



fj(6,u,t) = -^2f(9 + 2kTr/q,u},t)exp(2mjk9/q), 



fe=0 



(38) 



for j = 0, . . . , q. These modes satisfy the symmetry relation 
f j (0 + 2*/q,w,t)=exp(2nij0/q)f j (e,w,t). 

In analogy with the q = 2 state, we will find that the mode 
j = 0, corresponding to the symmetric part of / when q = 
2, has a globally-attracting low-dimensional description that 
evolves independently from the other modes, leaving q — 1 
free parameters to describe the distribution. 



A. Dynamics of the j = Mode 

A similar variation of the OA ansatz can be used to find 
a low-dimensional description of the dynamics of the j = 

I 



mode dynamics. The ansatz 



W,oo,t) = 



2tt 



l + J2a n (u},t)e qme + c.c.) , (39) 



yields the following ODE for a: 

a + q f iua + — (r q a 2 - r*) j = 0. 



(40) 



As before, we let g(u>) be Lorentzian with zero mean and 
spread A, such that r q (t) — a*(—iA,t) = a*(t), which 
closes the dynamics for r q 



|r |e**«: 



-A|r,| + y|r,|(l-|r,| a ) 



% = 0. 



(42) 



Thus, the manifold for the j : = mode dynamics, which can 
be shown to be globally attracting 0201 12111 . is the set of q-tuple 
Poisson kernels P(q8 — tp q ,\r q (t)\,ui). Again, we assume 
without loss of generality that %jj q = 0. 

B. Steady-State Solution 



With q potential clusters, the order-parameter \r q \ measures 
the degree of cluster synchrony in the system, while the lower 
q — 1 order parameters |n|, . . . , \r q -i\ measure the degree of 
asynchrony. Note that the distribution is only perfectly sym- 
metric if rx = • • ■ = r q _i = 0. Thus, there are q — 1 different 
measures of the asymmetry. 

Using a similar analysis as in the q = 2 case, we find that 
at steady-state 



rwM 



9(u) E.-=o(l/3 + CjW - <t>{") - 2kir/q) 



,)$? 



A-|r,|8in(^)]| 



if\oj\<K\r q 
if M > K\r a 



3in (*) /«• 



with \r q \ = y/1 - 2A/K and (j>(u) = 
Note that for f ss to be a distribution the coefficients Cj must 
satisfy Cj > —1/q and X^=o c j = 0' leaving q — 1 free pa- 
rameters that define the distribution. Note that in the q = 2 
case there was a single parameter [i.e., c in Eq. ( f20b l that char- 
acterized the asymmetry between the two clusters. 



The steady-state order parameters can be calculated using 
the same methods that led to Eq. d25l l. and analogous expres- 
sions (not presented here) can be obtained. 



r 



C. Transient Dynamics 



(43) 



To capture the transient dynamics, we study the PDE and 
corresponding characteristics given by 

d t f + [w-K\r g \ sin(q9)}def = qK\r q \ cos(q9)f, (44) 

=> = w-.fi:|j-g|siii(g0), (45) 

f = qK\r q \cos(q9), (46) 



-A\r q \ + ^\r q \(l-\r q f) 



(47) 



When \r q \ is at steady-state, we can solve Eqs. (|43T > and (l46l > 
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time t 

FIG. 9: (Color online) Comparison of \n\, |r"2 1, and [7*3 1 from the 
PDE J44t (blue circles, green triangles and red squares) to simulation 
with N — 10000 oscillators (cyan and magenta dashed lines) with 
the same initial conditions and parameters from figure[6] 



exactly, yielding equations analogous to Eqs. dAU and (IA2b in 
Appendix A for the characteristics of 6 and solution /, which 
we do not present here. 

When \r q \ is not at steady state its evolution is given by 



|r,(0l = P,l 



\ 



Pa 



e q(2A-K)t 



(48) 



where P q = y/l — 2A/K and Eq. (f4~4b can be solved nu- 
merically. In Fig. [9] we compare |n|, \r 2 \, and \r%\ from the 
numerically-computed PDE solution (blue circles, green tri- 
angles, and red squares, respectively) to a numerical simula- 
tion of Eq. d36*b with q = 3 and N = 10000 oscillators (cyan, 
yellow, and magenta dashed lines, respectively). The analytic 
solution for [t"3 | is plotted as a dot-dashed black line. 



V. DISCUSSION 

We have found an analytic description of both steady-state 
and transient dynamics of a system that shows cluster syn- 
chrony given by Eqs. (Q3 and (0. In the large N limit, q = 2 
solutions can be decomposed into symmetric and antisymmet- 
ric parts. The symmetric part, which evolves independently 
from the antisymmetric part and toward a steady-state inde- 
pendent of initial conditions, can be found using a variation 
of the OA ansatz 11911 and is globally attracting. The antisym- 
metric part, however, is shaped by the evolution of the sym- 
mertic part, is strongly dependent on initial conditions, and 



has a large degree of multistability. 

We have demonstrated how to manipulate the degree of 
asymmetry in the cluster states through the application of a 
short duration forcing term and modulation of the coupling 
strength. Starting from a symmetric state, any degree of asym- 
metry can be established by choosing the appropriate duration 
and strength of the forcing term. Furthermore, reducing the 
coupling strength decreases the amount of asymmetry in the 
cluster configuration, while increasing it does not have the op- 
posite effect, as shown in Fig. [HJa). Therefore, modulations 
of the coupling strength can be used to "erase" information. 
While we demonstrated this procedure using a system with 
q = 2, similar methods could be employed for q > 2. In 
particular, q — 1 parameters describe the cluster configura- 
tion, and the system could be driven to a configuration that 
encodes desired values of these parameters by the application 
of appropriately chosen forcing functions. Using these tech- 
niques, it is possible to encode information in the state of the 
system, which might find applications in the development of 
Kuramoto-type neural models. 

Problems that remain open include generalization such as 
the presence of noise and coupling functions with two or more 
harmonics. Thus far the work of Ott and Antonsen IU9I1 has 
not been generalized to these cases and no low dimensional 
analytic solution has been found. However, we hypothesize 
that when noise is added to Eq. ([3]) spontaneous switching can 
occur. The case where the coupling function has more than 
one harmonic has also been considered fiUl . In certain cases, 
e.g. H nm (9) = (K 1 sm(q 1 9) + K 2 sm(q 2 e))/NwhersK2 > 
K\, the resulting system is well-approximated to the class of 
systems studied in this paper and results, such as clustering 
and asymmetry, are qualitatively similar. 
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Appendix A: Characteristics 

In this appendix we present the results of solving the PDE 
in Eq. ( f2Tb via the method of characteristics when |r*2] is at 
steady-state (i.e. \r 2 \ = \/l — 2A/K). The characteristic 
ODEs are Eqs. d28l i and J29] i. Given an initial phase do, the 9 
characteristics evolve as 



J 



0(t,8o) = arctan 



K\r 2 \ — \/oj 2 — K 2 \r 2 \ 2 tan I arctan 



K\rn I —a; tan 9q 
^/fj 2 -A' 2 |r 2 | 2 



ty/u 



K 2 \r 2 \ 



LC 



(Al) 



Several example characteristics for the locked and drifting 
populations (w = 1 and cj = 3, respectively), are plotted in 



Fig.g|(a)and(b). 
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For initial conditions f(8, u>, £rj) = g(cu)h(6) the 9 charac- 
teristics can be used to solve for f(9, uj, t), given by 



m u,,t)=g(u;)-^-B D , 



(A2) 



where B 



K\r 2 \ sm[29(t)], and 



D= uj 2 + K 2 \r 2 \ 2 cos (E)- K\r 2 \\/u r 



K 2 \r 2 \ 2 sin (E) 



-K\r 2 \sin[26(t)} 



jK 2 



''2 



(A3) 



where E = 2 arctan 



(K\r 2 \ -ujttme{t))/y/Z? 
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